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Abstract. The problem of optimal mass transport arises in numerous ap- 
plications including image registration, mesh generation, reflector design, and 
astrophysics. One approach to solving this problem is via the Monge-Ampere 
equation. While recent years have seen much work in the development of 
numerical methods for solving this equation, very little has been done on the 
implementation of the transport boundary condition. In this paper, we propose 
a method for solving the transport problem by iteratively solving a Monge- 
Ampere equation with Neumann boundary conditions. To enable mappings 
between variable densities, we extend an earlier discretization of the equa- 
tion to allow for right-hand sides that depend on gradients of the solution 
[Froese and Oberman, SI AM J. Numer. Anal, 49 (2011) 1692-1714], This 
discretization provably converges to the viscosity solution. The resulting sys- 
tem is solved efficiently with Newton's method. We provide several challeng- 
ing computational examples that demonstrate the effectiveness and efficiency 
(O(M) - 0(M 1 - 3 ) time) of the proposed method. 



1. Introduction 

In this article, we propose a method for solving the elliptic Monge-Ampere equa- 
tion for a convex function u on a domain X £ M. d subject to the transport condition 
Vu : X — >• Y where Y is a connected set in M d . The method involves solving a 
sequence of Monge-Ampere equations with Neumann boundary conditions. We 
also describe an efficient finite difference method for solving these sub-problems. 
To demonstrate the capabilities of this solution method, we present computational 
results for several challenging numerical examples, which include the recovery of 
inverse maps, mapping onto unbounded density functions, mapping from a discon- 
nected domain, and mapping onto non-convex sets. 

1.1. L 2 optimal transport. The motivation for this work is the problem of opti- 
mal mass transport Amb03, Eva99 ( Vil03]. The problem originally considered by 
Monge is how to transport a given pile of sand into a hole in the most cost efficient 
way. Monge originally considered a cost equal to the magnitude of the distance the 
sand must be transported. More generally, the Monge-Kantorovich mass transport 
problem is to find a mapping s(x) that takes the density f(x) in the space X € M. d 
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to the density g(y) in the space Y £ K d and that minimizes the cost functional 

I[s]=fc( X ,s( X )) dx 

X 

where c(x, y) denotes the cost of transporting a unit of mass from the point x £ X 
to the point y £ Y. Here the data must satisfy the condition that total mass is 
conserved: 

J f(x)dx = J g(y)dy. 

X Y 

The problem of optimal mass transport arises in a number of important appli- 
cations including; image registration jHTKOll IHZTA041 luRHP+09) . mesh genera- 
tion jDCF+081 IFDC08I IBWOQj . reflector de sign [GOO gl IC()04| . and astrophysics 
(estimating the shape of the early universe) [FMMS02] . 

Kantorovich contributed to the understanding of optimal transport by refor- 
mulating the problem as a linear program and describing a simple dual formula- 
tion [Kan42, Kan48 . While this has made many theoretical questions easier to 
answer, this approach also effectively doubles the dimension of the problem. Con- 
sequently, computing the solution to even a small-scale problem is prohibitively 
expensive. This motivates the development of more sophisticated methods that 
will enable the efficient computation of optimal maps. 

In the special case of the quadratic cost function 

c(x,y) = - \x- y\ 2 , 

the problem has a special structure. In this situation, the optimal mapping s(x) can 
be expressed as the gradient of a convex function Eva99, Roc66 . The problem of 
obtaining the optimal mapping is then equivalent to the problem of solving a fully 
nonlinear partial differential equation (PDE) known as the elliptic Monge- Ampere 
equation 

(MA) det(D 2 u(x)) = f(x)/g{Vu(x)), x£X 

subject to the transport condition 

(BC) V« : X -> Y 

and the convexity constraint 

(C) u is convex. 

Remark 1 . Any solution method for the Monge- Ampere equation must enforce the 
convexity constraint, which is necessary to ensure a unique solution. 

1.2. Related works. In the past few years, the numerical solution of the Mongc- 
Ampere equation has received quite a bit of attention. However, most of the avail- 
able methods enforce Dirichlet, Neumann, or periodic boundary conditions rather 
than the transport condition that arises in many applications. 

An early work by Oliker and Prussner [OP88 presented a method that converges 
to the Aleksandrov solution of the Monge- Ampere equation in two dimensions. An- 
other convergent two-dimensional method was described by Oberman [O be08] ; this 
discretization converges to the viscosity solution of the equation. Other recent 
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methods, which perform best when solutions are sufficiently regular, have been de- 
veloped by Dean and Glowinski |DG061 lDG08l IGlo09j and Feng and Neilan |FN09al 
IFNOQbj . 

Recently, the author, together with co-authors, has extended the work of Ober- 
man |Obe06[ IObe08j to construct finite difference solvers that converge to the vis- 
cosity solution in any spatial dimension lil OK) F Olla) IFOf lb] . These methods 
perform quickly even in the most singular examples. 

Much less work has been done on the implementation of the transport bound- 



ary condition (BC). A fluid flow approach was introduced by Benamou and Brc- 
nier [BBOOj and has been further developed by Haber, Rehman, and Tannen- 
baum [HRT10 . However, this approach is computationally expensive as it requires 
introducing an additional dimension to the problem. We also mention the work by 
Finn, Delzanno, and Chacon |FDC08j , which enables the mapping of a square to a 
region with four (possibly curved) sides. For the related problem of optimal trans- 
port (or partial transport) with cost given by the distance c(cc, y) = \x — y|, a finite 
element method has been constructed by Barrett and Prigozhin |BP07|, [BP09 . 

1.3. Contents. In |scction 2| of this work, we review some analysis — including weak 
solutions and regularity results for the Monge- Ampere equation — that inform the 
approach taken in this paper. In |scction 3| we describe our method for implement- 



ing the transport boundary conditions. In section 4 we provide a discretization 



of the Monge-Ampere equation and prove that it converges to the viscosity solu- 
tion. In |scction 5} we provide further details about the numerical implementation 
of our method. In|scction 6| we provide computational results that test our Monge- 



Ampere solver. In sectioii 7} we provide computational results for several challeng- 



ing and representative transport problems. In section 8 we summarize the main 
contributions of this work. 



2. Analysis and weak solutions 

In this section, we review some regularity results for the Monge-Ampere equation 
that are needed to fully explain the approach taken in this work. 

2.1. Weak solutions. The Monge-Ampere equation is a second order PDE, so 
classical solutions of this equation should have at least two continuous derivatives. 
However, these classical C 2 solutions do not always exist. Thus it is necessary to 
use some notion of weak solution to make sense of non-smooth solutions of the 
Monge-Ampere equation. 

Weak solutions of the Monge-Ampere equation can be defined in different ways. 
The discretization used in this work is motivated by the viscosity solution (which, 
in most cases, is equivalent to the more general Aleksandrov solution). 

We recall the definition of viscosity solutions [CIL92j , which are defined for the 
Monge-Ampere equation in [G ut01| . 

Definition 1. Let u G C(X) be convex and f > be continuous. The function 
u is a viscosity subsolution (supersolution) of the Monge-Ampere equation in X if 
whenever convex <p € C 2 (X) and xq € X are such that (u — <j))(x) < (>)(it— (f))(xo) 
for all x in a neighborhood of Xq, then we must have 



detp 2 ^)) > (<)f(x ). 
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The function u is a viscosity solution if it is both a viscosity subsolution and su- 
persolution. 

Example 1 (Viscosity solution of Monge- Ampere) . For concreteness, we provide 
a particular example of a function that, though not a classical C 2 solution of the 



Monge- Ampere equation, can be understood as a viscosity solution. Consider (MA) 
with solution and / given by 

U (x) = i((|x|-1)+) 2 , /(x) = (l-l/|x|)+. 



This function u (pictured in Figure 1 ) is a viscosity solution — but not a classical 
C 2 solution — of the Monge- Ampere equation. 

We verify the definition of a viscosity solution. This only needs to be done at 
points where |x | = 1 (since u is locally C 2 away from this circle). We note that / 
is equal to zero on this circle. 

We begin by checking convex C 2 functions <fi < u with <^>(xo) = it(xo) = (that 
is, u — <p has a local minimum here). Since Vu(xo) = 0, we require V0(xo) = 
as well. Since u is constant in part of any neighborhood of Xo, any convex 4> must 
also be constant in this part of the neighborhood in order to ensure that u — cf> 
has a local minimum. This means that </> has zero curvature in some directions, so 
that detD 2 cj)(xo) = 0, as required by the definition of the viscosity solution. We 
conclude that u is a supersolution of the Monge- Ampere equation. 

We also need to check functions <fi > u with </>(xo) = u(xo) = (so that u — <fi 
has a local maximum) . Since <p is convex, it will automatically satisfy the condition 
det D 2 <fi(x ) > and we conclude that u is a subsolution. 




2.2. Regularity. We now review the regularity we can expect for solutions of the 
L 2 optimal transport problem. These results are due to Caffarelli |Caf92al [Caf92b , 
ICafflB] . 

We begin by noting that, as with the Monge- Ampere equation, solutions of the 
transport problem need not be smooth. An example of a singular solution (see 



Figure 2 1 is the problem of mapping the circle 



X = {(£1,0:2) I x\ +x\ < 1} 
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onto the disconnected set 

Y = {(x 1 ,x 2 ) | xi < -0.25, (x-l + 0.25) 2 + x\ < 1} 

U {(xi,x 2 ) | X! > 0.25, (x t - 0.25) 2 + x\ < 1}. 

In fact, the solution remains singular even if the disconnected region Y is ap- 
proximated by a connected region Y e . 

While we do not solve the problem of mapping onto a disconnected region, we 
are able to solve for the inverse mapping (which takes the disconnected set Y to 



the connected set X) in £ 7.2 




As long as the sets X,Y are bounded, we are at least guaranteed that the solution 
of the Monge-Ampere equation is differentiable almost everywhere with bounded 
gradient. 

Remark 2. When the solution to the Monge-Ampere equation is not differentiable, 
the map is given by the sub-gradient rather than the gradient. This allows a single 
point to be mapped onto a region rather than a single point. 

More regularity is guaranteed if we restrict ourselves to convex target sets Y. 

Theorem 1 (Interior Regularity). Suppose thatX,Y are bounded, connected, open 
sets and Y is convex. Suppose also that the density functions 

f:X-> (0,+oo), g : Y -> (0, +oc) 

are bounded away from and +oo . Then the solution of the Monge-Ampere equa- 
tion {MA}, jBC] ), |[C|| belongs to C]£{X) for some < a < 1. 

//, in addition, the density functions f, g G C' 9 for some < /3 < 1 then the 
solution of Monge-Ampere belongs to Cf^(X) for every < a < (3. 

If both sets X, Y are uniformly convex, we can obtain regularity up to the bound- 
ary as well. 

Theorem 2 (Boundary Regularity). Suppose, in addition to the hypotheses of 
\Theorem 1\ that the sets X and Y are uniformly convex. Then the solution of 
Monge-Ampere is in C 2 ' a {X) for some < a < 1. 
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3. Transport boundary conditions 

In this section, we discuss the transport boundary conditions in more detail. We 
describe a method for solving this challenging problem by solving a sequence of more 
tractable sub-problems; these are Monge-Ampere equations subject to Neumann 
boundary conditions. 



3.1. Nonlinear boundary conditions. In the problem of L 2 optimal transport 
between convex sets X, Y <E R rf , the transport condition (BC) 

Vw : X -> Y, 

also known as the second boundary value problem, can be enforced by simply 
requiring boundary points to map to boundary points |Pog71[ ITW091 IUrb97j : 

Vu : dX -> BY. 

In particular, if the boundary of the region Y is defined by the function 

<%) = 0, 

we can write the transport boundary condition as 
(1) $(V«(x)) =0, x e dX. 

While we might try simply enforcing this nonlinear equation at boundary points, 
the function <f> can be highly nonlinear and non-smooth. As a result, it will be 
difficult to construct a discretization that is consistent with a possibly singular 
solution of the equation and that will permit fast solvers to remain stable. 



3.2. Mapping to rectangles. The situation simplifies significantly if we are sim- 
ply mapping a rectangle to a rectangle. In this case, since the optimal L 2 mapping 
does not permit twisting or rotation, we expect the four sides of the rectangle X 
to map to the corresponding si des of the rectangle Y . 

As a concrete example (see Figure 3), suppose that the sets X, Y € R 2 are 
defined as 

X = (0,1) x (0,1), Y= (0,1) x (0,1). 

Then, for example, we expect the function Vu(x) to map the segment x\ = 0, x 2 € 
[0, 1] to the segment yi — 0, y 2 £ [0, 1]. That is, 

«xi (0,2:2) = 0. 

Similarly, we will have 

u Xl (l,x 2 ) = 1, u X2 (xi,0) = 0, u X2 (x 2 , 1) = 1. 



This is simply a (linear) Neumann boundary condition, which is straightforward to 
implement [BS911 lQbe06] . 
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Figure 3. Mapping between squares. 



3.3. A sequence of Neumann boundary conditions. Given the appearance 
of the gradient in the transport boundary condition ([I]) and the simplicity of im- 
plementing a Neumann boundary condition, we would like to find the Neumann 
boundary condition 

^ = 4>{x), xedx 

on 

for the Monge- Ampere equation that is equivalent to solving the more challenging 
problem (MA), (BC), |C]). Here the vector n refers to the unit outward normal 



vector at each point x € dX. 

It is not at all apparent from (BC I what the equivalent Neumann boundary con- 
dition should be. However, we suggest a sequence of Neumann boundary conditions 
that can be used to numerically determine the correct function <fi. 

We first recall that the gradient of the exact solution u maps the boundary of 
the set X to the boundary of Y 

Vu:dX^ dY 

and that the correct Neumann condition is given by 

4>(x) = Vu(x) ■ n(x), x € dX. 

To find this function, we suppose that we have a convex approximation u k to the 
solution of the Monge- Ampere transport problem. Then the (sub-)gradient of this 
function will map the domain X onto some set Y k G R d and, since u k is convex, 

Vu k : dX -> dY k . 

In reality, we would like the image of the gradient to be dY, the boundary of the 
target set. This motivates us to consider the projection of dY k = Vu k (dX) onto 
the correct set of boundary points dY: 

Yio\ 9Y {\/u k {x)) = argmin||y- Vw fe (x)||2, x e dX. 

yedY 

From this we extract a new Neumann boundary condition 

k (x) =Proj ar (Vu fe (a;)) -n 

and solve the Monge-Ampere equation once again with this updated boundary 
condition to obtain a new approximation u k+1 . 
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To summarize, we iterate to produce a sequence of functions (it 1 , 
by solving the Monge- Ampere equation 

'det(D 2 u k+1 {x)) = f(x)/g(Vu k + 1 (x)), 
X?u k+1 (x) ■ n(x) = Proj ar (Vu fc (x)) ■ n : 



.) obtained 



(2) 



4> k {x), 



xex 
xedX 



u k+1 is convex. 



We make the important observation that these boundary conditions do not pin 
down the values of Vw fc+1 on the boundary. This would be a mistake since we know 
only that Vm : dX — > dY and not the exact values of Vu(x) on the boundary. 
Instead, each Neumann condition fixes only one component of the gradient (the 
normal component) and allows the remaining component(s) to slide as needed to 
ensure that the Monge- Ampere equation is satisfied. 

3.4. Solvability of sub-problems. We note that the iteration ^ may not be 
well-posed. The problem here is that, while the Monge- Ampere equation with the 
correct Neumann values <j)(x) has a solution, the sub-problems we have described 
may not be solvable. 

We recall that for the Monge- Ampere equation with a Neumann condition: 

'det{D 2 u) = f{x)/g{Vu{x)), xeX 
Vu(x) ■ n(x) — ip(x), x € dX 



u is convex, 



x e X 
xedx 



a solution (unique up to an additive constant) exists only if an implicit solvability 
condition is satisfied |LTU86) . 

To get around this problem, we instead solve a problem of the form 

(det(D 2 u) = cf(x)/g(Vu(x)), 
\ Vu(x) ■ n(x) = ip(x), 
I u is convex, 
(J x u(x) dx = 

for the unknowns c > and u[x), where the constant c is chosen to ensure the 
equation has a solution and the mean-zero condition forces the solution to be unique 
(instead of unique up to an additive constant). 

Of course, if we are given the correct Neumann values 4>(x) for the solution to the 
transport problem, the constant c will simply be equal to one. However, by relaxing 
this condition we make it possible to solve the sub-problems when the solvability 
condition requires c to be slightly different than one. 

To summarize, we solve the transport problem by performing the iteration 

' det{D 2 u k+1 {x)) = c k+1 f{x)/g{Vu k+1 (x)) 1 x G X 

X?u k+1 (x) ■ n(x) = Proj aF (Vu fc (2 ; )) 1 n i x ) = 4> k { x ), x G dX 



(3) 



u k+1 is convex, 
J x u k+1 (x) dx = 0. 



Remark 3. As we pointed out in {2.2 solutions to the optimal transport problem 



need not be continuously differentiable up to the boundary. In this case, the Neu- 
mann boundary condition should be understood in the viscosity sense [Lio85J. To 
obtain the boundary condition (f> k (x) at a point where the gradient Vu k {x) is not 
defined, we instead look at the projection of a value in the sub-gradient of u k {x). 
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4. DISCRETIZATION 

Having described an iteration for solving the transport problem, we now need 
to describe the method for solving the sub-problems. The most challenging step 
here is to discretize the Monge- Ampere equation. In this section, we describe a 
discretization that provably converges to the viscosity solution. 

4.1. Standard finite difference methods. The simplest thing to do is to simply 
discretize the Monge- Ampere equation using standard centered differences. 

In two dimensions, for example, the Monge- Ampere equation has the form 

Ux 1 x 1 Ux 2 x 2 ~ u x\x 2 = f( x )/9( u x 1 ,U X2 )- 

A standard centered difference discretization of this equation is 

(4) MA h s [u] = (V XlXl u){Vx 2 x 2 u) - {V XlX2 u) 2 - f{x)/g{V Xl u,V X2 u) 

where the finite difference operators are defined by 

[V X2X2 u]ij = — (ui >j+1 + Uij-i - 2uij) 

[Dx 1 x 2 u ]ij — T7^2 + l + u i-l,J-l ~ + l — u i+lj-l) 

[V Xl u]i] = 2h _ u i-l,j) 

\Dx 2 u\ l3 = — (U ilj+1 - UiJ-l) . 

However, as is pointed out in FOllb , this discretization may fail to converge 
to the correct viscosity solution and solution methods can become unstable. 

4.2. Convergent finite difference methods. Before we discuss our discretiza- 
tion of the Monge-Ampere equation, we briefly review the theory of convergent 
finite difference methods for viscosity solutions of nonlinear elliptic equations. The 
foundation for this is a result by Barles and Souganidis |BS91| . 

Theorem 3 (Convergence of Approximation Schemes). Consider a degenerate el- 
liptic equation, for which there exist unique viscosity solutions. A consistent, stable 
approximation scheme converges uniformly on compact subsets to the viscosity so- 
lution, provided it is monotone. 

Oberman |Obe06| used this result to further characterize convergent finite dif- 
ference discretizations. We recall that a finite difference equation has the form 

F l [u] = F l (u,,u, - Uj\jjti). 

Then a degenerate elliptic (monotone) scheme can be defined as follows: 

Definition 2. The scheme F is degenerate elliptic if F l is non- decreasing in each 
variable. 

Theorem 4 (Convergence of Finite Difference Discretizations). Consider a degen- 
erate elliptic equation, for which there exist unique viscosity solutions. The solution 
to a consistent, degenerate elliptic finite difference scheme converges uniformly on 
compact subsets to the viscosity solution. 
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(a) In the interior. (b) Near the boundary. 

Figure 4. Wide stencils on a two dimensional grid. 



Even for linear elliptic equations, it is not always possible to construct monotone 
schemes using a narrow stencil, as was demonstrated in the early work by Motzkin 
and Wasow [MW53j . Oberman |Obe081 used wide stencils to construct monotone 
discretizations of second directional derivatives for directions v lying on the grid. 
These derivatives can be discretized using centered differences: 

(5) T> uu Ui = — R — { u ( x i + vh) + u { x i ~ vh) — 2u(xi)) . 

\v\ h 2 

Depending on the direction v, this may involve a wide stencil. At points near 
the boundary of the domain, some values required by the wide stencil will not 
be available ( Figure 4"| . In these cases, we use interpolation at the boundary to 



construct a (lower accuracy) stencil for the second directional derivative; see |Obe08j 
for more details. The consistency error of these approximations depends on both 
the spatial resolution h and angular resolution d9 of the stencil. 

4.3. Discretization of the Monge-Ampere operator. Next we describe a mono- 
tone discretization of the Monge-Ampere operator 

det(£> 2 u), 

which was introduced in |F011aj . This discretization is a consequence of the fol- 
lowing variational characterization of the Monge-Ampere operator 

d 

det(D 2 u) — min I I max{w„.„., 0} 

where u is a convex function and V is the set of all orthonormal bases of R d : 
V = {fa, ...,v d ) | Vj G M. d , Ui -L Vj if i £ J, \Wj\\ 2 = !}• 

We are going to add an additional term to this expression in order to further 
penalize non- convexity: 

!d d 
TT max{it„.„. , 0} + \] ^^{UuAUi , 0} 
i=i 

We note that if u is a convex function, the terms in the summation will vanish since 
all second directional derivatives u v . v . are non- negative. On the other hand, if u 
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is a non-convex function then at least one of the directional derivatives u vv will be 
negative and we will have 

det+(D 2 w) < u uv < 0. 
Thus this non-convex function cannot be a solution of the Monge- Ampere equation 



(6) det + (L» 2 u) = F(x, Vu) > 0. 

In order to discretize this, we restrict our attention to a finite set of orthogonal 
vectors Q on the grid. A monotone discretization of the Monge- Ampere operator 
is given by 

!d d 
TT max{Z>„.„.'u ) 0} + V" vahx{V VjV ,u, 0} 
3 = 1 

where T> vv u is the discretization of the second directional derivative that is defined 
in ^. 

4.4. Discretization of functions of the gradient. In the works [FOlla , FOllb], 
we considered problems where the right-hand side / was a function of x only. 
However, in the Monge-Ampere equation that arises in the transport problem, 
the right-hand side also depends on the gradient of the solution. That is, the 



P Xj u(x) = — (w(x + hej) — u(x — hej)) 



equation ( MA I is of the form 

(7) det(D 2 u(x)) = F(x,\7u(x)). 

Consequently, we need to discretize not only the eigenvalues of the Hessian but also 
the gradient. 

The simplest approach would be to simply use standard centered differences for 
the first derivatives: 

1 

2h' 

where e^ is the vector whose i th component is equal to the Kronecker delta 8ij. 
While this discretization is consistent with C 2 solutions of the Monge-Ampere 
equation, it is not monotone and there is no guarantee that it will converge to 
the viscosity solution. 

Oberman |Ob e06 provided some examples illustrating the construction of mono- 
tone discretizations for functions of the gradient. For example, that work describes 
a monotone discretization of the absolute value of a first derivative: 
1 
h 

For more general functions of the gradient, one approach to producing a mono- 
tone discretization is to simply use centered differences and add on a small multiple 
of the laplacian: 

g(u x ) = g{V x u) + hK g V xx u + 0{h). 

Here K g is the Lipschitz constant of the function g. 

However, instead of adding an additional term to the discretized equation, we 
could make use of the second derivatives that are already present in the Monge- 
Ampere equation. This is the subject of the following section. 



I^K^i)] = t max{u(i£j) — u(xj-i), u{xj+\) — u{xj), 0} + 0{h). 
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4.5. Discretization of the Monge-Ampere equation. So far we have attempted 
to produce a monotone discretization for each individual term in the Monge-Ampere 
equation. As an alternative to this, we suggest using a wide stencil to produce a 
discretization of the Monge-Ampere equation which, though it may not be mono- 
tone for each of the individual terms, is monotone when considered as a whole. This 
discretization also ensures that the linear systems that must be solved in the im- 
plementation of Newton's method are well-conditioned even when the eigenvalues 
of the Hessian are close to zero or the iteration is initialized poorly. 

To accomplish this, we make use of the second directional derivatives u Vj Vj that 



are already present in the Monge-Ampere equation, as noted in \ 4.4 By making a 
change of coordinates, we can write the gradient 

Vw = (u Xl ,. . .,u Xd ) 

in terms of first derivatives in the directions vy. 

V« = {u Vl ,. ..,«„„). 

Once this is done, the only problem we might have will be if one of the second 
derivatives vanishes, in which case the Monge-Ampere operator will not be uni- 
formly elliptic. We can remedy this by simply regularizing the maximum and 
minimum functions slightly to bound them away from zero: 

max{-, 0}, min{-, 0} — > max{-, <5}, min{-, 5} 

where S > is a small parameter. 

To accomplish all this, we first need to rewrite the gradient in terms of the new 
coordinate system. We consider any set of d orthogonal vectors in M d : (vi, . . . , vj). 
Now we can rewrite the gradient of a function u in terms of directional derivatives 
along these axes: 



d d 



Vu= (u Xl ,...,u Xd ) = \J2^u7T Uv ^---'Yl 



,Vi\ — \ v i 

3=1 ' 3 ' 3=1 ' J 



This enables us to discretize the gradient using a wide stencil by discretizing the 



directional derivative in the direction Vj as 



1 



(8) V V] Ui = ^ (u(xi + vjh) - u(xi - Vjh)) , 

which has an accuracy of 0{h 2 ). Near the boundary, where some of the required 
values may not be available, we can simply use a first-order accurate forward or 
backward difference. We stress again that this discretization of the gradient is valid 
for any set of orthogonal vectors v%, . . . , vj- 
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Using this characterization of the gradient, we can rewrite the Monge- Ampere 
equation as 

!d d "I 

TT max{u v . v . , 0} + min {u v . v . , 0} > - F(x, Vu) 
7=i U J 

!d d 
TT waxUivtv ,, 0} + y^min {u v . v .,Q \ - F(x, Vu) 

!d d 
TT vobx{u v . v . , 0} + min {u u , v , , 0} 

'ft N ft N 

= (, 1 , m S) 6 v G (-- 

As we have already described in ([5|,([8]), the directional first and second deriva- 
tives can be discretized using a wide stencil by limiting the set of possible directions 
in the set V to a finite set Q of orthogonal vectors that lie on the grid. We also 
introduce a small parameter 8 > in order to bound the maximum and minimum 
functions away from zero: 

max{-, 0}, min{-, 0} — > max{-, 6}, min{-, 5}. 

We can now define the discretization of the Monge- Ampere equation as 

(9) MA h 1 f> s [u] = t min_G^>] 



(»i,..,» d )eg to.-."*) 1 

,,1 io AaG 

(ui,...,U4) I 



where each of the Q^' de ' S \ u ] [ s defined as 



G (^f.'Sd) M = II max-fZ^.u, 5} + ^ mm{V VjUj u, 5}- 

W / d 

I ft N ft N 

Theorem 5 (Convergence to Viscosity Solution). Let the PDE ^ have a unique 
viscosity solution and let the right-hand side F(x, Vu) be Lipschitz continuous on 
fl x R d with Lipschitz constant Kp. Then the solutions of the scheme ^ converges 
to the viscosity solution of (MA| as h,d8,6 — >• with S^ 1 > Kp \vj\ h/2 for every 

vj e g. 

Proof. The convergence follows from verifying consistency and degenerate elliptic- 
ity. This is accomplished in Lemmas [7]|8] □ 



Lemma 6. Under the hypotheses of Theorem 5, the scheme for Q^ de > s [ u ] { n ( io ) 
is degenerate elliptic. 

Proof. We introduce the notation 

pt(xi) = u(xi + huj) - u(xi), pj(xi) = u(xi - huj) - u(xi). 
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This allows us to write Q^d ,S \ u ] [ n the form of Definition 2 as follows: 

(V 1 ....,v d ) L J |_| 

(ii) GbZt«M^> ■ ■ ■ >pt,p~ d ) = n max (f^ff 

+ Vnfin{^^,4-^f^^^--4 :: ff 

Now we need only check that this is non-decreasing in each of its arguments. We 
verify this for the term p^; the reasoning is identical for the remaining terms. 
Choose any e > and consider: 

Gh,d9,S i — I — i \ /-th.dO.S / 4-\ 

>^-i Lax { V ± + \ + Pi A - max / g£+gL , ; V) 

V I N ^ 2 J In* 2 jy 

+ L n M +e , + ft ,4 - min { £ + s\) 



hi 2 /* 2 1 I kil 2 ^ 2 



P i + e -Pi _ Pi ~Pi 
2\v 1 \h ~ 2\v 1 \h 

In the above, we have used the facts that 



min J pt+i+K A min J pt+K A > o 

I W ^ 2 J lM 2 * a J" 

and that <5 < 1. 

We continue with this expression to conclude that 

rd-i f Pi" + e +?r , r Pi" +pr x 



> S a - L 11 • ^ + S - ^ Y -5\- K F —— 
{S*- 1 - K F \ Vl \h/2). 



This expression is positive as long as 6 > Kp \ v\\h/2. 

We conclude that each of the G^' d9 ' > is increasing in each of its arguments 
and is thus degenerate elliptic. □ 

Lemma 7. Under the hypotheses of Theorem^ the scheme for MA^jf 6 ' 6 ^] in |9]) 
is degenerate elliptic. 

Proof. This scheme is the minimum of degenerate elliptic schemes, and is therefore 
degenerate elliptic. □ 

Lemma 8. The scheme for Al A^f 6 ' S [u] in ([9| is consistent with the Monge- Ampere 
equation (|6| for any function u G C 2 (X) . 

Proof. The proof of this is identical to the consistency proof in |FQllai Lemma 6] . 

□ 
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4.6. Hybrid discretization. As in |F011b| . we can improve the accuracy of the 
discretization by using the monotone scheme only in regions of the domain where 
the solution may be singular. In smooth regions of the domain, we simply use a 
standard centered difference discretization. 



Given the regularity results described in {2.2 we cannot expect solutions to be 
smooth at the boundary since the domain X that we are computing on is a square, 
which is not strictly convex. However, we can expect more regularity in the interior 
of the domain as long as the density functions are C a and are bounded away from 
and oo. 

We first identify X s , which is a neighborhood of any regions where u may be 
singular: 

X s = dXU{x £ X | f(x) < e}L){x £ X \ f(x) > l/e}L){x eX\f$ C a {B(x,e))}. 

Here e is a small parameter, which we can take equal to the spatial resolution h. 

We define w(x) to be a function that is one in an h- neighborhood of X s and that 
goes to zero elsewhere. Then a possible hybrid discretization of the Monge- Ampere 
equation is 

(12) MA h ' M ' 5 [u] = w{x)MA h jf ,5 [u] + (1 - w{x))MA h s [u]. 

5. Numerical implementation 

Now that we have described the discretization we will be using, we turn our 
attention to the remaining details of the numerical implementation of the iteration 
described in this paper. In this section, we describe our methods for enforcing 
boundary conditions, solving the discrete system of equations, and initializing the 
iteration. To be concrete, we will describe these issues in the two-dimensional case, 
but the ideas easily generalize to higher dimensions. 

5.1. Existence and uniqueness. We recall here that the iteration ^ requires us 
to solve not only for the function u, but also for the scaling factor c that multiplies 
the density functions. We can simply include this as an additional variable in the 
discrete system of equations. 

With the addition of an extra variable, we should expect that we will also require 
an additional equation. This is reasonable since solutions of the Neumann problem 
are only unique up to an additive constant. We simply add an extra equation that 
forces u to be zero at one corner of the domain. 

5.2. Boundary conditions. We also need to discretize the Neumann boundary 
conditions. 

Our computational domain is the square, which means we must impose values 
for u Xl on the left and right sides of the domain and for u X2 on the top and bottom 
edges of the domain. 

We accomplish this by adding a layer of ghost points around the outside of our 
computational domain. The value of the normal derivatives on the boundary can 
then be discretized using simple centered differences. For example, at a point on 
the left edge (x% = x m - m ), we can discretize the normal derivative as 

U a (x) = ^-(w(x min + h,x 2 ) - u(x min - h,x 2 )). 

The use of ghost points ensures that all values needed in this discretization are 
available. 
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We also need to provide four more equations at the corner points in our grid. We 
specify the value of the derivative in the "diagonal" direction ((1,1), (1,-1), (—1,1), 
or (— 1, —1)) that points outward from the grid at each of these four points. This is 
enforced using centered differences. So, for example, at the points (xi ; mm, ^2,min) 
we require that 

(u(a;i,min - h, X 2 ,min ~ h) - (Xl, m in + K ^2,mm + k)) = 



2y/2h 

^2 (^^1 (^l,min, ^2,min) ~l~ (^1 , mill 3 *£2,min) ) ■ 

As before, the ghost points ensure that all of these values are available. 

5.3. Newton's method. The discretization of the Monge-Ampere equation de- 
scribed in Sj4] results in a system of equations that can be solved efficiently using 
Newton's method. 

This involves performing the iteration 

u k+1 = u k - v k , c k+l =c k - d k 

where the correctors v k , d k are obtained by solving the equation 

\7MA[u k , c k }(v k , d k ) T = MA[u k , <*]. 

As long as the initial iterate u° satisfies the given Neumann boundary condition, 
we can simply enforce a homogeneous Neumann condition on the corrector v k at 
each step. 

Although we are using a hybrid discretization, the weight function that deter- 
mines the discretization is independent of the iterates u k , c k . This means that 
we can compute the Jacobians of the monotone and standard discretizations and 
obtain the Jacobian of the hybrid system via 

VMA[u,c] = w(x)VMA M [u,c] + (1 - w(x))VMA s [u, c]. 

We begin by computing the Jacobian of the monotone discretization. We recall 
that this discretization has the form 

MA M [u,c]= min G(u 1 ,...,u d )[u, c]. 
(ui,...,v d )eg 

By Danskin's Theorem [B er03j . we can write the Jacobian of this as 

VM4 M [ u ,c] = VG (l , lr . w) [M,c], 

where the (y\, . . . , is<i) are the directions active in the minimum. 

This Jacobian can be broken down into two basic components: the gradient with 
respect to the solution vector u and the gradient with respect to the scaling factor 
c. The first component is given by: 



V Ui G (l/1) ... )!/(i) [ti,c] = 
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The final component is given by 

d d 

v c G il/ly .. tVd) [u,c] = -f ( x ,Y J \^ v u i n,---,Yy\^r LVv i Ui 




We also require the Jacobian of the standard discretization. In this case, the 
first component of the Jacobian (in two dimensions) is simply 



V Ui MA s [u,c] = (V 

x2x 2 u i)T^x 1 x 1 + (T^x 1 x 1 Ui)T) X2X2 + 2(T) XlX2 iii)T) XlX 

dF dF 

c^—(x,T> Xl u i ,T) X2 u i )'D Xl - c-—(x,T> Xl u i ,T> X2 u i )'D X2 
dpi - dp 2 

and the second component is 



V c Mi s [u,c] = -F(x,T> Xl Ui,V X2 Ui). 

5.4. Initialization. The iterations we have described in this paper also need to 
be initialized. There are really two aspects to this: we need to initialize u and c 
each time we solve the Monge- Ampere equation and we also need to initialize our 
estimation of the boundary conditions <j>(x). 

5.4.1. Initialization of boundary data. First we discuss the initialization of the 
boundary data cf>° in the iteration ([3]). The simplest approach would be to extract 
boundary conditions from the identity map s(x) = x. However, if this mapping 
does not overlap with the target set Y, the iteration is likely to fail. 

We can remedy this problem by instead extracting boundary data from the scaled 
identity map s(x) = Mx where the constant M is chosen large enough so that the 
set s(X) encompasses the target set Y. 

Once this constant is chosen, we simply choose the initial boundary condition 

(j)°(x) = Mx-n(x), xedX. 

We can accelerate the convergence of this method by first solving the transport 
problem on a coarser grid, then interpolating the resulting boundary data onto the 
refined mesh. 



5.4.2. Initialization of Newton's method. We also need to initialize Newton's method 
each time we solve the Monge- Ampere equation. We can use the approach described 
in [FOllaj . which involves obtaining the initial guess by solving the equation 

Au{x) = (cd\f(x)/g(x-x ))V d 

where xq is a point in the interior of the target set Y. 

Since we will be solving the Monge- Ampere equation multiple times with differ- 
ent boundary conditions, we can also accelerate the convergence of the (k + l) st 
iteration by initializing with the solution found during the previous solve (u k ). One 
important point here is that the boundary data changes from step to step. Thus 
it is important to change the values of u k at the boundary points so as to ensure 
that correct boundary conditions are satisfied. 
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6. Computational results: the Monge-Ampere equation 

In this section, we provide computational results for several different examples. 
We have really introduced two ideas in this paper: a discretization for Monge- 
Ampere type equations and a method for implementing the transport boundary 
condition. In this section, we focus on testing our Monge-Ampere solver. To 
keep this idea clear, we restrict ourselves to the problem of mapping rectangles 
to rectangles; in this case, our method for the transport problem is reduced to a 



single Monge-Ampere solve with Neumann boundary conditions (see f 3.3 ) . 

In each example our domain is a square, which is discretized on an N x N grid 
using a 17 point stencil. We let h — 1/(N — 1) denote the spatial resolution of the 
grid and let M = N 2 denote the total number of grid points. The computations 
were done in MATLAB on a laptop with a 2 GHz Intel processor. 

When an exact solution u exact is available, we provide the maximum error in the 
gradient map: 

Error - max{||<f ct - u Xl \U ||<f ct - ^ 2 |U}. 

We also provide the total number of Newton iterations and computation time re- 
quired for each example. 

The examples we consider include: 

• A (linear) map between gaussian densities. 

• A comparison between a map obtained by solving the direct problem and 
a map obtained by inverting the solution to the inverse problem. 

• A map from a uniform density onto a density that blows up at a point. 

• A map between two brain MRI images. 



6.1. Gaussian densities. We begin by showing that we can recover a linear map- 
ping between two rectangles with gaussian densities. We consider the problem of 
mapping the square (—0.5,0.5) x (—0.5,0.5) onto the rectangle (0.5, 1.5) x (—1, 1) 
with the density functions: 



f(xi,x 2 ) = ^ ex P 



r 2 
•'1 



1 

2 0.4 2 



2 0.4 2 



l{VUVl) 



1 

O08 



exp 



1 (m i) 2 i 

2 0.4 2 2 0.2 2 



In this case, we have an explicit expression for the optimal map: 

1 



u Xl = xi + 1, 



-x 2 . 



We present the results in |Tablc"T and Figure 5 In this example, we can actually 
achieve machine accuracy (if we take enough Newton steps). This is because the 
exact solution is simply a linear map, which will exactly solve the discretized system 
of equations. In addition to this, we find that the Newton solver for the Monge- 
Ampere equation converges in O(M) time. 
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N h Newton Iterations CPU Time (s) Maximum Error 



10" s 
10" 8 
1(T 8 
10" 8 
10~ 8 
10" 8 
10" 8 
10" 8 

Table 1. Computation time and maximum error for the map be- 



32 
46 
64 
90 
128 
182 
256 
362 



0.0323 
0.0222 
0.0159 
0.0112 
0.0079 
0.0055 
0.0039 
0.0028 



0.1 
0.2 
0.3 
0.6 
1.1 
2.4 
5.3 
12.4 



5.71 x 
3.34 x 
0.26 x 
0.18 x 
0.13 x 
0.09 x 
0.07 x 
0.05 x 



tween two gaussian densities (E 6.1 ) 



6.2. Recovering an inverse map. For our next example, we consider another 
problem with an exact solution, which will be used to verify that we can correctly 
recover inverse maps. To set up this example, we define the function 

q(z) = (-i-z 2 + ^ + cos(8ttz) + ^zsin(8^). 

Now we map the density 

f(x 1 ,x 2 ) = l+4(g''( Sl ) g (x 2 )+g( a ; 1 )( Z "( a ; 2 ))+16( g (x 1 )g(a ;2 ) (7 ''(a; 1 )< Z ''( a ; 2 )- ( 7'(x 1 ) 2 q'(x 2 ) 2 ) 

in the square (—0.5,0.5) x (—0.5,0.5) onto a uniform density in the same square. 
This transport problem has the exact solution 

u Xl {xi,x 2 ) = xi +4q'(xi)q(x 2 ), u X2 (xi,x 2 ) = x 2 + Aq(xi)q'(x 2 ). 

We will solve this problem in two ways: 

• Directly, as in the previous example. 

• By solving the inverse problem (mapping g to /) and inverting the resulting 
map. 

Results are presented in |Figurc 6] and |Table 2| We find that the maps obtained 
from both the forward and inverse formulations have about 0(h 2 ) accuracy. Both 
problems are solved in about O(M) time. 
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Forward Problem Inverse Problem 



N 


Iterations 


Time (s) 


Max Error 


Iterations 


Time (s) 


Max Error 


32 


3 


0.2 


2.476 x 10 -3 


4 


0.4 


2.450 x 10 3 


46 


2 


0.2 


0.631 x 10~ 3 


2 


0.5 


0.575 x 10 3 


64 


2 


0.5 


0.241 x 10~ 3 


2 


1.1 


0.244 x 10 3 


90 


1 


0.6 


0.106 x 10~ 3 


1 


1.3 


0.101 x 10 3 


128 


1 


1.3 


0.049 x 10~ 3 


1 


2.9 


0.048 x 10 3 


182 


1 


2.9 


0.024 x 10~ 3 


1 


5.1 


0.023 x 10 3 


256 


1 


6.3 


0.012 x 10~ 3 


1 


10.9 


0.011 x 10 3 


362 


1 


14.0 


0.006 x 10~ 3 


1 


22.6 


0.006 x 10 3 



Table 2. Newton iterations, computation time and maximum er- 
ror for a map obtained by a direct solve and by inverting the inverse 
map (Q. 



6.3. An example with blow-up. Next we consider the problem of mapping a 
uniform density onto a density that blows up at a point: 

exp (-2v/(yi - 0.5) 2 + (y 2 - 0.5) 2 ) 

0(2/1,2/2) = j. n = . n = • 

v/(2/i - 0.7) 2 + (y 2 - 0.7) 2 

In this case, both X and Y are the square (0, 1) x (0, 1). This example is taken 
from [DCF+08] . which allows us to compare results. In this example, we slightly 
regularize the density g (bounding it by a 0(l/h?) function) to prevent infinities 
from appearing. 

We present the timing results in Table [3] We provide not only the number of 
Newton iterations and computation time, but also the ratio 

R = m&x{g(y 1 ,y 2 )/ f(x 1 ,x 2 )} , 

since many currently available Monge- Ampere solvers can become slow or unstable 
when this ratio is large. For comparison, we provide the same information for the 
method of [DCF+08] (which is essentially our "standard" discretization solved with 
an optimized Newton-Krylov method). The method of DCF+08] runs in 0{M) 
time. Our method, though it runs in about 0{M ) time, has lower computation 
times and deals with larger density ratios. Naturally, we cannot conclude too much 
from the comparison of computation times since the computations were performed 
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on different computers. However, it is evident that, in terms of computation time, 
our method is very competitive with other fast solvers. 

We also present the deformed mesh and zoom into the region of high density to 
verify that our method has produced an untangled mesh; see Figure [7J 



Hybrid Method Method of [DCF+08 



N R 


Iterations 


CPU Time (s) 


R 


Iterations 


CPU Time (s) 


32 546 


4 


0.2 


356 


6 


1 


46 1,151 


4 


0.3 








64 2,254 


5 


0.8 


1,127 


7 


4 


90 4,066 


5 


1.6 








128 9,162 


5 


3.5 


2,829 


7 


17.4 


182 18,608 


5 


8.3 








256 36,933 


5 


19.4 


8,886 


7 


70 


362 74,018 


4 


36.3 








Table 3. 


Ratio of density functions, Newton iterations 


, and total 



computation time for the hybrid method (j |4.6|5.3 ) and the method 
of jDCF+08] , 




6.4. Mapping between brain MRI images. We conclude this section with an 
example from image processing. In this example, we obtain our density functions 
from the pixel intensities in two synthetic brain MRI images |Cenlf)l ICZK+98I 
CKKSE97J. The images are shown in Figures 8(a)|8(b) In this case, the regions 



X and Y are identical and are equal to the unit square. The fully resolved images 
contain 256 x 256 pixels. For the computations presented here, we have also in- 
terpolated both images onto coarser grids so that in each case we are mapping an 
N x N grid onto the density function obtained from an N x N image. 

In this example, the density functions have large gradients, which effectively 
increase as we map onto more refined images. The solver now runs in about 0(M 1A ) 
time; see |Table~4} 
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Figures |8(c)|8(d)| show the image we obtain by solving the Monge- Ampere equa- 
tion and interpolating and the error in this image. The mapped image we obtain 
agrees well with the given image. Not surprisingly, the largest error occurs around 
the edges of the brain where the density function is essentially discontinuous; conse- 
quently, small errors in the map can lead to large errors in estimated pixel intensity. 




Figure 8. 
function g, 



(c) (d) 
(a) | The initial density function /, |(b)| the final density 
(c)| the image obtained by solving the Monge- Ampere 



equation and interpolating, and|(d)|the error in the resulting image 



N Newton Iterations CPU Time (s) 



32 
4(3 
64 
90 
128 
182 
256 



7 
7 
9 

10 
12 
12 
13 



1.1 
1.2 

3.0 
7.0 
13.7 
34.9 
81.6 



Table 4. Computation time for a map between two brain MRI 
images ([6.4). 
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7. Computational results: optimal transport 

In this section, we turn our attention to computational results for the trans- 
port problem. In each example, we embed our domain in the square (—0.5, 0.5) x 
(—0.5,0.5) (setting the density / = outside our domain X). While this can lead 
to singularities in the solutions, our methods are robust enough to handle this 
non-smoothness. 

In each case, we present the total number of Monge-Ampere solves required on 
the N x N grid (this does not include solves performed on coarser grids during the 
initialization process), as well as the total computation time required. When an 
exact solution is available for comparison, we provide the maximum error in the 
map: 

Error = max{||<f ct - u^, \\u™ act - u^}. 
The examples considered in this section include: 

• A map between two ellipses, for which an exact solution is available for 
comparison. 

• A map from two disconnected semi-circles onto a circle, for which an exact 
solution is available for comparison. 

• A map from a square onto a convex polygon, which is neither smooth nor 
strictly convex, together with recovery of the inverse map. 

• A map from a square onto a non-convex region. 

7.1. Mapping an ellipse to an ellipse. First we consider the problem of mapping 
an ellipse onto an ellipse. To describe the ellipses, we let M x ,M y be symmetric 
positive definite matrices and let B\ be the unit ball in M. d . Now we take X = M X B\, 
Y = M y B 2 to be ellipses with constant densities /, g in each ellipse. 
In M. 2 , the optimal map can be obtained explicitly MOO fj from 

Vu(x) = M y R g M~ 1 x 

where Rg is the rotation matrix 

cos(0) -sin(0) 



Re ~ 1 sin(0) cos((9) 
the angle 9 is given by 

tan(6>) = tr&ceiM- 1 My 1 J)/trace(M^ 1 My 1 ), 
and the matrix J is equal to 

J = Rn/2 = ( 1 q 

We use the particular example 

0.4 \ / 0.3 0.1 



Mx ~ \ 0.2 ) ' My ~ \ 0.2 0.4 

which is pictured in |Figure 9| 

Projections onto the ellipse at each step are accomplished efficiently using the 
method described in |Kis94| . 

Computational results are presented in |Table 5] and |Figure 9] The error is 
decreasing uniformly (about O(h o s )). We cannot expect high accuracy for this 
example due to its degeneracy: the density / vanishes in part of the domain. This 
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means that the lower accuracy monotone stencil is needed in this region, which will 
in turn affect the error in the map. 

Despite the degeneracy of this example and the multiple Monge- Ampere solves 
required to initialize and solve this problem, the computation requires only 0(M 11 ) 
time. 




N h \MA\ Solves CPU Time (s) Maximum Error 

32 0.0323 4 0.7 0.0264 

46 0.0222 13 1.7 0.0180 

64 0.0159 3 1.8 0.0152 

90 0.0112 6 5.5 0.0117 

128 0.0079 3 9.9 0.0083 

182 0.0055 3 25.3 0.0060 

256 0.0039 2 61.9 0.0048 
Table 5. Computation time and maximum error for the map be- 
tween two ellipses (j|7.1[). 



7.2. Mapping from a disconnected region. We now return to the degenerate 
example considered in §2.2| This is the problem of mapping the two half-circles 

X = {(x!,x 2 ) | x x < -0.1, Oi +0.1) 2 +x 2 2 < 0.3 2 } 

U {(xi,x 2 ) | xi > 0.1, (x x - 0.1) 2 + x\ < 0.3 2 } 

onto the circle 

Y = {(yi, y 2 ) \ y\ + yl < 0-3 2 }. 

Results are presented in |Table 6| and |Figurc 10} In this case, the error appears 
to approach a constant value of around 0.004. This is not surprising since the 
monotone scheme is needed in the region where / vanishes or is discontinuous. The 
width of the stencil then limits the accuracy of solutions; this point is explained 
more fully in |FQllaj . The computation time for this very degenerate example is 
about C(M 13 ). 
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-0.2 

x 1 

(a) 



0.2 




(b) 



Figure 10. (a) A cartesian mesh in two half-circles X and |(b)| its 
image under the gradient map Vw ([7.2). 



N 


h 


( MA 1 Solves 


CPU Time (s) 


Maximum Error 


32 


0.0323 


5 


0.5 


0.0171 


46 


0.0222 


2 


0.5 


0.0160 


64 


0.0159 


5 


1.6 


0.0129 


90 


0.0112 


9 


6.0 


0.0082 


128 


0.0079 


5 


11.8 


0.0052 


182 


0.0055 


4 


30.3 


0.0040 


256 


0.0039 


3 


66.7 


0.0038 



Table 6. Computation time and maximum error for the map from 



two half-circles to a circle ([ 7.2 ) 



7.3. Mapping to a convex polygon. Next we consider a map onto a convex 
polygon Y , which has a very non-smooth boundary. We use the polygon Y with 
vertices: 

(-0.5,-0.3), (-0.5,0.4), (0,0.5), (0.5,0.3), (0.3,-0.5). 

Despite the non-smoothness of dY, our method successfully maps the square (—0.5, 0.5) x 
(—0.5,0.5) into the prescribed polygon, though we do not have an exact solution 
to compare with. 

We also compute this map by solving the inverse problem (mapping the polygon 



to the square) and inverting the map as in [ 6.2 While no exact solution is available 
for comparison, we can check the maximum difference between components of the 
two maps: 

rii tnv ii ii mv w ~\ 

max{\\u Xl - u Xl Woo, \\u X2 - u X2 
Results are presented in |Tablc 7 and Figure 11 The computation is reasonably 



efficient, requiring about 0(M 12 ) time for both the forward and inverse problem. 
We also observe that the agreement between the maps obtained from the forward 
and inverse approaches improves as we refine the grid. 

7.4. Mapping to a non-convex region. Finally, we compute the mapping of the 
square with constant density / onto a non-convex region given by 

Y = {(yi,y 2 ) | < Vl < 1, < y 2 < 1 - 0.1sin(27r yi )} . 
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0.5 



-0.5 




(a) 




(b) 



Figure 11. (a) A cartesian mesh and |(b)| its image under the 
gradient map Vw (j]7.3[). 



N 


Forward Problem 
Iterations Time (s) 


Inverse Problem 
Iterations Time (s) 


Max Difference 


32 


3 


0.4 


1 


0.3 


0.0397 


46 


3 


0.8 


1 


0.7 


0.0227 


64 


3 


1.5 


1 


1.1 


0.0153 


90 


4 


3.2 


1 


2.3 


0.0119 


128 


4 


8.5 


1 


6.2 


0.0087 


182 


4 


21.0 


1 


13.5 


0.0063 


256 


4 


61.8 


1 


33.9 


0.0050 


362 


4 


154.3 


1 


92.6 


0.0044 



Table 7. Monge- Ampere solves, computation time and maximum 
difference for a map from square to polygon obtained by a direct 



solve and by inverting the inverse map ({ 7.3 1 



We impose the following periodic density in the region Y: 

<7(yi,J/ 2 ) = 2 + cos (SttV^j - 0.5) 2 + (y 2 - 0.5) 2 ) 



The results are displayed in Table 8| and |Figurc 12| Despite the non-convexity of 
Y, the method successfully maps the region X into the non-convex region Y . The 
non-convexity does not appear to affect the computation time at all: the solution 
time is roughly O(M). 



N I MA I Solves CPU Time (s) 



32 
46 
64 
90 
128 
182 
256 
362 



1.7 
2.2 
5.4 
8.4 
21.4 
41.9 
68.1 
197.4 



Table 8. Computation time for the map onto a non-convex region ([7.4 1 
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8. Conclusions 

In this paper, we proposed a numerical method for solving the Monge-Ampere 
equation with transport boundary conditions. The method we described requires 
the solution of a sequence of Monge-Ampere equations with Neumann boundary 
conditions. In practice, we observed that for sufficiently large problems, the number 
of Monge-Ampere solves was independent of or even decreased with the size of 
the problem. No more than 13 iterations were required in any of the examples 
considered in this work. 

We also presented a discretization of the Monge-Ampere equation that extends 
the work of [FOllaj to the more general setting where the right-hand side can de- 
pend on the gradient of the solution. We proved that the solution of the discretized 
system converges to the viscosity solution of the equation and described Newton's 
method for solving the resulting system efficiently 

Finally, we provided computational results for several challenging and represen- 
tative examples including recovering inverse maps, mapping onto an unbounded 
density, mapping onto a non-convex target set, mapping from a disconnected do- 
main, and mapping between images. In each case, our method successfully and 
efficiently (O(M)-O(M 13 ) time) computed a map into the specified target set. 
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